# bcp.R
# R code for Bayesian change point estimation described in the paper,
# Bruce Western and Meredith Kleykamp, "A Bayesian Change Point Model for
# Historical Time Series Analysis." The paper describes a Gibbs sampler
# and several approaches to stochastic sampling from the conditional
# posterior (SSCP) of regression coefficients.
# Bruce Western, April 27, 2004

library(MASS)

#--------------------------------------------------------------
# Defining some useful functions
#--------------------------------------------------------------
# Linear predictor function
lp <- function(x,b) as.vector(apply(t(x)*b,2,sum))

# Normal likelihood function
n.like <- function(eta,y,sigma) {
	ll <- sum(dnorm(y,mean=eta,sigma,log=TRUE))
	exp(ll) }

# Determinant of a matrix
det <- function(x) prod(eigen(x)$values)

#-------------------------------------------------------------
# Gibbs function
#-------------------------------------------------------------
gibbs <- function(X, y, tstar=length(y)/2, iter=500, burn=100) {
# User supplies X, a matrix including unit vector in col. 1 for intercept
#               y, a dep var vector
#               tstar, start value for the change-point
#               iter, number of iteration in the gibbs sampler
#               burn, number of burn-in iterations

# OLS fit
out <- lm(y~X-1)
p <- ncol(X)
n <- length(y)

# Setting up data
grid1 <- (p+1):(n-p-1)

XX <- as.list(grid1)
for(i in grid1) {
	tt <- (1:n)>=i
	XX[[i-min(grid1)+1]] <- cbind(X, X * tt)
}

# Priors
pb <- rep(0,2*p)
pvb <- diag(rep(1e6,2*p))
pvbi <- diag(1/rep(1e6,2*p))
ps2 <- .001
pdf <- .001

# Start values
bstar <- rep(out$coef,2)

# Initializing parameter vectors
beta <- h <- tt <- NULL
k <- length(grid1)
temp <- rep(NA,k)
thetas <- grid1

# Posterior variance quantities
posta <- (n/2)+pdf

for(i in 1:(iter+burn)) {
	indt <- as.numeric((1:n)>=tstar)
	XZ <- cbind(X, X*indt)
	xx <- crossprod(XZ)
	xy <- crossprod(XZ,y)
	yhat <- lp(XZ, bstar) 
	estar <- y-yhat
	posts2 <- ps2+sum(estar^2)            
	hstar <- rgamma(1,posta,rate=posts2/2)
	postvb <- solve(pvbi+(xx*hstar))
	postb <- as.vector(postvb %*% (xy*hstar))
	bstar <- as.vector(mvrnorm(1,postb,postvb))
	temp <- unlist(lapply(lapply(XX,lp,bstar),n.like,y,sqrt(1/hstar)))
	postt <- temp/sum(temp)
	tstar <- sample(thetas, 1, prob=postt)
	if(i>burn){
		beta <- rbind(beta,bstar)
		h <- c(h,hstar)
		tt <- c(tt,tstar)
} }
mode <- table(tt)
mode <- names(mode)[mode==max(mode)]
output <- list(grid1,tt,h,beta,as.numeric(mode))
names(output) <- c("time","theta","precision","beta","mode of theta")
output
}

#-------------------------------------------------------------
# SSCP function
#-------------------------------------------------------------
sscp <- function(X, y, pmu, pvb, pa, pb, bayes=T) {
# User supplies: X, a matrix including unit vector in col. 1 for intercept
#                y, a dep var vector
#                pmu, vector of prior means
#                pvb, prior covariance matrix for coefs
#                pa, prior degrees of freedom
#                pb, prior sums of squares
#                bayes=T, coefs based on bayes posterior probs
#                bayes=F, coefs based on llik posterior probs

n <- length(y)
p <- ncol(X)

# Grid over which break point is searched
ind <- (p+1):(n-p-1)

tau <- solve(pvb)*pb/pa

betam <- betapv <- ppr <- as.list(ind)
ppm <- ll <- rep(NA,length(ind))
postb <- NULL
astar <- pa + (n/2)

for(i in ind) {
	indt <- as.numeric((1:n) >= i)
	xm <- cbind(X, indt*X)
	ols <- lm(y ~ xm - 1)
	yhat <- predict(ols)
	sigma <- summary(ols)$sigma
	ll[i-p] <- n.like(yhat, y, sigma)
	pxpx <- crossprod(xm) + tau
	betam[[i-p]] <- solve(pxpx) %*% (tau%*%pmu + crossprod(xm, y))
	e <- as.vector(y - (xm%*%betam[[i-p]]))
	b2 <- as.vector(crossprod(pmu - betam[[i-p]],tau%*%pmu))
	dm <- pb + sum(e*y)/2 +  b2/2
	ppr[[i-p]] <- (astar/dm) * pxpx
	ppm[i-p] <- dm^(-astar) * det(pxpx)^(-.5)
}
ppm <- ppm/sum(ppm)
ppll <- ll/sum(ll)
output <- cbind(ppm,ppll)
dimnames(output) <- list(as.character(ind),c("Bayes","LL"))

# SSCP for coefficients
if(bayes) {
	mixind <-  table(sample(ind,10000,replace=TRUE,prob=ppm))
	mode <- ind[ppm==max(ppm)] }      
else  {
	mixind <- table(sample(ind,10000,replace=TRUE,prob=ppll))
	mode <- ind[ppll==max(ppll)] }      
newind <- as.numeric(names(mixind))
for(i in 1:length(newind)) {
	mu <- as.vector(betam[[newind[i]-p]])
	Sigma <- solve(ppr[[newind[i]-p]])
	normb <- mvrnorm(mixind[i], mu, Sigma)
	tb <- normb/sqrt(rchisq(1,astar)/astar)  
	postb <- rbind(postb,tb) } 

output <- list(ind, output, postb, as.numeric(mode))
names(output) <- c("time","p(theta|y)","beta","mode of theta")
output
}

